Homework 4

Name:

Collaborators:

Due Friday. October 30th 2015

Question 1 : Solving Linear Systems

In this homework you will implement one of the most basic linear solvers. Gaussian elimination with bakwards substitution with row swap, you can find the details in Alg 6.1 in your textbook (or https://en.wikipedia.org/wiki/Gaussian_elimination)

Gaussian Elimination

Q1.a You will write the Gaussian elimination algorithm. We will proceed by writing the function rowSwap!. The inputs of the function will be:

  • A, a square matrix
  • i,j the indeces of the rows to be swaped.

The function will not have any output. The swaping has to be done in the matrix M. Finally, your function will raise an error if $i$ or $j$ are out of range.

Remark:

In Julia, all the input variariables are passed by reference. Then you can modify the input variables outside the scope of the function. By convention, each time that a function modifies a variable outside its scope the function contains an exclamation mark.


In [7]:
function rowSwap!(A, i,j)
    # input:   A matrix
    #          i,j the row indexes to be swapped
    n = size(A)[1]
    # checking that i and j are within the permitted range
    (i > n || j > n ) && error("Index out of range in the rowSwap")
    # if the index are not the same
    if i != j 
        buffer = A[i,:]
        A[i,:] = A[j,:]
        A[j,:] = buffer
    end
end


Out[7]:
rowSwap! (generic function with 1 method)

Q1.b You will write a function that performs the Gaussian elimination. The function takes as input:

  • A a square matrix,
  • b a vector.

Your function will create the augmented system, and perform the Gaussian elimination. The output of the function will be the tuple (U,b1). U is the upper triangular matrix resulting from the elimination and b1, is the resulting vector.
To obtain $U$ and $b1$ from your augmented matrix you perform a slicing (i.e. use [:,1:n]).
You will use your rowSwap function acting on the augmented matrix.
Finally, your function will raise an error if:

  • the matrix is not square,
  • the matrix is singular,
  • the dimension of the vector and matrix are not compatible.

Hints:

  • You may use the function 'hcat' to build the augmented matrix. In Julia you can type
    ? hcat
    and press enter to obtain some information of the function.
  • To find the pivot you can use the function "find()".
  • If the vector has only zero elements then the function "find()" will output an empty vector.
  • You can check that your matrix is non-invertible by checking that the output of find is empy.
  • You can easily check if an vector is empty using the "isempty()" function with a true or false output.

In [21]:
function gaussianElimination(A,b)
    #input:   A squared matrix
    #         b a vector
    #output:  U upper triangular matrix
    #         b1 the resulting vector 
    
    # safety checks
    (n,m) = size(A)
    (n != m) && error("Matrix is not square \n")
    (n != size(b)[1]) && error("Dimension mismatch \n")
    
    # create the augmented matrix 
    M = hcat(A,b)
    for i = 1:(n-1)
        # look for the first non zero entry in the ith column
        # find the indices of the non zero elements (up to machine precision)
        indexj =  find(abs(M[i:end,i]).> eps(1.0))
        # if indexj is empty then the matrix is singular and raise error
        isempty(indexj) && error("The matrix is singular \n")
        # call row swap
        rowSwap!(M,i,(i-1)+indexj[1])
        # for loop for eliminating unknows
        for j = i+1:n
            M[j,:] = M[j,:] - (M[j,i]/M[i,i])*M[i,:]
        end
    end
    # checking the last pivot!! 
    abs(M[n,n]) < eps(1.0) && error("The matrix is singular \n")
    # slicing the matrices
    U = M[:,1:n]
    b1 = M[:,n+1:end]
    return (U,b1)
end


Out[21]:
gaussianElimination (generic function with 1 method)

Triangular solver (backwards substitution)

Once the matrix is reduced to a upper triangular form, the system can be solved by backsubstitution.

Q1.c You will write a function that performs the backward substitution.
The input of your function will be:

  • U: an upper triangular matrix,
  • b: a vector.

The output will be the solution to the system $Ux = b$.
Your function needs to have safeguards against a size mismatch (i.e., the size of the matrix and your vector are not compatible, or your matrix is not squared).

Hint: if you need to run a foor loop that goes from n-1 to 1, you can use the syntax
for j = n-1:-1:1


In [9]:
function backwardSubstitution(U,b)
    # input:    U upper triangular matrix 
    #           b vector
    # output:   x = U\b
    # checks for sizes
    (n,m) = size(U)
    (n != m) && error("Upper triangular matrix is not square \n")
    (n != size(b)[1]) && error("Dimension mismatch \n")
    
    # creating x and running the backward substitution
    x = zeros(b)
    x[n] = b[n]/U[n,n]
    for i = (n-1):-1:1
        x[i] = (b[i] - dot(U[i,i+1:end][:],x[i+1:end]))/U[i,i]
    end
    
    # returning x
    return x
end


Out[9]:
backwardSubstitution (generic function with 1 method)

You can test that your function is correct by running the following script:


In [10]:
# size of the Matrix
m = 100
# creating an upper triangular Matrix 
U = diagm(m*ones(m,1)[:], 0) + diagm(rand(m-1,1)[:], 1) + diagm(rand(m-2,1)[:], 2)
# creating the rhs
b = rand(size(U)[1],1)
@time x = backwardSubstitution(U,b)
print("Residual of the backward substitution is ", norm(U*x -b)/norm(b),"\n")


  

The residual should be extremely small (around epsilon machine)


In [11]:
function backwardSubstitutionLoop(U,b)
    # checks for sizes of U and b
    n = size(U)[1]
    # first step
    x = zeros(b)
    x[n] = b[n]/U[n,n]
    for i = (n-1):-1:1
        dummy = 0
        for j = i+1:n
            dummy += U[i,j]*x[j]
        end
        x[i] = (b[i] - dummy)/U[i,i]
    end
    return x
end


Out[11]:
backwardSubstitutionLoop (generic function with 1 method)

Linear solver

Q1.d You will write a function (very short) that solves a linear system in the form $Ax = b$, using your function GaussianEliminiaton and backwardSubstitution.
The input of your function will be :

  • A, a square matrix
  • b a vector.
    The output will be the answer $x$.

In [25]:
function solveGauss(A,b)
    # input:    A square matrix 
    #           b vector
    # output:   x = A\b
    (U,b1) = gaussianElimination(A,b)
    return backwardSubstitution(U,b1)
end


Out[25]:
solveGauss (generic function with 1 method)

You can test your code by running the following script


In [26]:
# size of the Matrix
m = 100
# creating the Matrix 
A = randn(m,m) + m*eye(m)
println("The conditioning number of A is " ,cond(A))
# creating the rhs
b = rand(size(A)[1],1)

@time x = solveGauss(A,b)
print("Residual of the solver is ", norm(A*x -b)/norm(b),"\n")


The conditioning number of A is 1.3136771554235003
  0.010945 seconds (39.43 k allocations: 18.396 MB)
Residual of the solver is 4.767478691948936e-16

In [30]:



The conditioning number of A is 1
Out[30]:
11x11 Array{Float64,2}:
 -0.342327     0.0108213     0.660125     …   0.134672     -0.268872   
  0.0          4.12741e-12   6.1483e-12      -1.72959e-12   1.3686e-12 
  0.0          0.0          -6.24985e-13     -1.11183e-12  -3.4754e-13 
  0.0          0.0           0.0              5.86496e-13  -1.5901e-13 
  5.55112e-17  0.0           0.0             -3.72681e-12   1.2657e-12 
 -8.31066e-17  0.0           0.0          …   8.34487e-13  -3.16339e-13
  2.77979e-16  0.0           0.0             -1.03499e-14  -1.27151e-12
 -1.51874e-16  0.0          -5.04871e-29      3.714e-12    -3.56918e-12
 -2.78615e-17  0.0           4.49572e-29     -1.2428e-12    9.38704e-13
 -1.18404e-16  0.0          -1.81762e-28      7.04903e-13  -1.0759e-12 
 -1.00659e-18  0.0          -8.72815e-29  …   5.04871e-29  -4.9665e-13 
.2390180313431834e18
  0.010552 seconds (36.25 k allocations: 18.234 MB, 25.84% gc time)
Residual of the solver is 1.4760558826374621

The residual should be extremely small (around epsilon machine)

Question 2: Complexity

From your textbook you know how many operations you need to solve the a linear system based on

You will perform a benchmark to obtain the assymptotic complexity of your algorithm with respect to the number of unknows in the system.
You will run your algorithm to solve linear systems for matrices of different sizes; you will time the execution time for each of those solves, and you will plot the runtime versus the size of the matrices in a log-log scale. From the plot you will claim the assymptotic complexity of your algorithm with respect to the number of unknowns on the linear

Q2.a You will run the following script, to bechmark the assymptotic complexity of the Julia built-in linear system solver (which is an interface to LAPACK, for more information see https://en.wikipedia.org/wiki/LAPACK and http://www.netlib.org/lapack/)


In [14]:
nSamples = 10;
times = zeros(nSamples,1)
sizes = 2*2.^(0:nSamples-1)
for i = 1:nSamples
    m = sizes[i]
    # creating the Matrix 
    A = rand(m,m) + m*eye(m)
    # creating the rhs
    b = rand(size(A)[1],1)
    tic(); 
    x =A\b
    times[i] = toc();
end


0.020930 seconds (39.39 k allocations: 18.398 MB, 11.90% gc time)
Residual of the solver is 6.104396625319271e-11
elapsed time: 1

You will plot the timing using the following script (you will use the layer object to plot two different data-set in the same graph)


In [15]:
using Gadfly
plot(
    layer(x = sizes, y = times, Geom.point,  Geom.line),
    layer(x = sizes, y = 0.00000001*sizes.^3, Geom.point, Geom.line, Theme(default_color=color("red"))), 
    Scale.y_log10, Scale.x_log10,
    Guide.ylabel("Runtime [s]"), # label for y-axis
    Guide.xlabel("Size of the Matrix"),  # label for x-axis
)


LoadError: ArgumentError: Gadfly not found in path
while loading In[15], in expression starting on line 1

 in require at ./loading.jl:233

In red you have the cubic scaling and in light blue the numerical runtimes. What should you expect?
What are you seeing instead?
How can you explain this?
What would happen if you increse the size of the matrices?

Answer: We should expect the two lines to be parallel; however, the zig-zag pattern look like a quadratic complexity. This is a clear example of a prea-assymptotic regime, given that the problem is not big enough we are not able to "see" the leading exponent. If the size of the problems to solve were bigger, we would be able to see the assymptotic scaling.

Q2.b You will modify the script in the question above in order to time the algorithm you wrote.


In [16]:
nSamples = 10;
times2 = zeros(nSamples,1)
sizes2 = 2*2.^(0:nSamples-1)
for i = 1:nSamples
    m = sizes2[i]
    # creating the Matrix 
    A = rand(m,m) + m*eye(m)
    # creating the rhs
    b = rand(size(A)[1],1)
    tic(); 
    x = solveGauss(A,b)
    times2[i] = toc();
end


.08028453 seconds
elapsed time: 2.9409e-5 seconds
elapsed time: 4.0305e-5 seconds
elapsed time: 5.7436e-5 seconds
elapsed time: 0.00339448 seconds
elapsed time: 0.000247837 seconds
elapsed time: 0.00073755 seconds
elapsed time: 0.001610172 seconds
elapsed time: 0.004690021 seconds
elapsed time: 0.02889815 seconds
elapsed time: 3.3517e-5 seconds
elapsed time: 2.6259e-5 seconds
elapsed time: 7.6071e-5 seconds
elapsed time: 0.000123385 seconds
elapsed time: 0.000552901 seconds
elapsed time: 0.00204418 seconds
elapsed time: 0.019955855 seconds
elapsed time: 0.120047632 seconds
elapsed time: 0.982154334 seconds
elapsed time: 14

In [17]:
using Gadfly
plot(
    layer(x = sizes2, y = times2, Geom.point,  Geom.line),
    layer(x = sizes2, y = 0.00000001*sizes2.^3, Geom.point, Geom.line, Theme(default_color=color("red"))), 
    Scale.y_log10, Scale.x_log10,
    Guide.ylabel("Runtime [s]"), # label for y-axis
    Guide.xlabel("Size of the Matrix"),  # label for x-axis
)


LoadError: ArgumentError: Gadfly not found in path
while loading In[17], in expression starting on line 1

 in require at ./loading.jl:233
.249635337 seconds

Based on the runtime scaling you obtained, what is the assymptotic complexity of your function solveGauss?

We can see the two lines are parallel (at least for the last 3 or 4 points) then we can claim that the complexity of the algorithm is cubic as expected.


In [ ]: